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1 Introduction 

1.1 Objective 

This report summarizes work performed by the Space Systems Laboratory (SSL) for NASA 
Langley Research Center in the field of performance optimization for systems subject to 
uncertainty. The objective of the research is to develop design methods and tools to the 
aerospace vehicle design process which take into account lifecycle uncertainties. It recognizes 
that uncertainty between the predictions of integrated models and data collected from the 
system in its operational environment is unavoidable. Given the presence of uncertainty, 
the goal of this work is to develop means of identifying critical sources of uncertainty, and 
to combine these with the analytical tools used with integrated modeling. In this manner, 
system uncertainty analysis becomes part of the design process, and can motivate redesign. 
The specific program objectives were: 

1. To incorporate uncertainty modeling, propagation and analysis into the integrated 
(controls, structures, payloads, disturbances, etc.) design process to derive the error 
bars associated with performance predictions. 

2. To apply modern optimization tools to guide in the expenditure of funds in a way that 
most cost-effectively improves the lifecycle productivity of the system by enhancing 
the sub-system reliability and redundancy. 

The results from the second program objective are described in [21]. This report de- 
scribes the work and results for the first objective: uncertainty modeling, propagation, and 
synthesis with integrated modeling. 

1 . 2 Report Overview 

This report is broken into four parts, covering the areas researched. Section 2 provides 
background on uncertainty in system design as investigated in the SSL. It includes a litera- 
ture search performed on propagation of uncertainty, and describes possible approaches to 
dealing with uncertain system interfaces. Section 3 uses knowledge on the sources of uncer- 
tainty along with data taken from the Multi-body Active Control Experiment (MACE) to 
experimentally derive uncertainty bounds for several parameters. Section 4 describes the 
change-of- variables propagation routine developed by SSL. Finally Section 5 integrates all 
of the research together and shows an example problem of how to use uncertainty analysis 
in the design process. 

This integration is performed in the DOCS analysis environment. DOCS (Disturbance- 
Optics-Controls-Structures toolbox), developed at MIT and commercialized by Mide Tech- 
nology Corporation is a MATLAB® toolbox developed for the creation and dynamical 
analysis of integrated models (figure 1). Integrated aerospace models can be created using 
finite element program, control systems, and disturbance or optics models. These models 
are first numerically conditioned, and then can be run through a host of suites. These 
include predicting the root-mean-square (RMS) values of performance metrics, sensitivity 
analyses which can be coupled with gradient-based optimization, and uncertainty analysis 
to include error bars with the predicted performance values. The research described by this 
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report has used and built upon DOCS; specifically the boxes ’’Uncertainty Database” and 
” Uncertainty Analysis” in figure 1 were advanced. 
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Figure 1: Roadmap for the DOCS framework. Red boxes indicate that new capabilities 
have been added. 


1.3 MACE 

The Multi-body Active Control Experiment (MACE) is an example of an aerospace system 
to illustrate the use of uncertainty in design. MACE (figure 2) is a 1.6 meter long flex- 
ible structure designed and built in cooperation by MIT, Payload Systems Incorporated. 
Lockheed Martin Missiles and Space Company, and NASA Langley Research Center. It 
was originally built to investigate high precision pointing and vibration control of flexible 
structures [19], and has been used in developing techniques for system identification, model 
updating, and nonlinear control. 

MACE consists of four struts connected axially to each other, with aluminum nodes at 
each interface. Three of the struts are polycarbonate, and the fourth is an active strut with 
piezoelectric actuators. At the center of this boom is a reaction wheel assembly consisting of 
three momentum wheels. At each end are gimbals in the X- and Z-directions, which provide 
disturbance inputs on one side, and control actuation on the other. Below the gimbals are 
0.7 meter long flexible appendages. 

MACE comes fully instrumented with a complete actuator /sensor package. There 20 
sensors comprised of rate gyros, optical encoders, strain gages and reaction wheel tachome- 
ters. The 9 actuators are made up of gimbals, three reaction wheels, and piezo-electric 
active struts. A separate Electronics Support Module (ESM) provides real time control, 
and any necessary actuator and sensor conditioning. 
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Figure 2: Illustration of the MACE testbed. 


2 Uncertainty in Systems Design 

The challenge in validating high performance systems is that modeled performance pre- 
dictions cannot be independently validated by systems level tests until they are in their 
operational environment, at which point it may be extremely difficult and costly to make 
gross changes to the design. 

An alternative to systems testing is to validate the design using modeling, analysis, and 
lower-level testing. An illustration of this approach is shown in Figure 3. Since uncertainty 
will continue to exist in all models, its presence must also be taken into consideration. As 
yet, uncertainty in system design is not treated in a rigorous fashion. Nominal values of 
model results are tracked, but results are frequently given without any confidence intervals. 
Uncertainty safety factors are often used in place of any real knowledge of the system 
uncertainties. This highlights a need in current design practices to include uncertainty in 
the design process. 

The goal of this research is to build a framework for uncertainty management in the 
design process for high performance systems. The research includes identifying model un- 
certainties, evaluating uncertainties to determine how much they affect the system perfor- 
mance, and tracking both sources of uncertainty, or inputs to the model, and overall system 
uncertainty, or model outputs, through the process. An overview of literature in the subject 
follows. 

2.1 Literature Overview 

There is a wealth of literature on the subject of uncertainty. Many papers have been writ- 
ten which use simple methods to place tolerance bounds about the nominal value for some 
metric. Researchers in such diverse fields as engineering, economics, and political science 
have applied uncertainty analysis techniques to models of everything from aerospace struc- 
tures, bridges, and chemical processes, to economic trends and socio-environmental conse- 
quences [8]. All of these works tend to use one of a small handful of uncertainty techniques. 
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Figure 3: Design lifecycle paradigm, as developed at the MIT Space Systems Laboratory. 


Below is a brief overview of the techniques, along with advantages and disadvantages. 

The great majority of uncertainty analyses in literature deal with the propagation of 
an uncertain input (or inputs) through a model in order to determine the uncertainty of 
an output. These individual model propagation techniques can be grouped into bounded 
or probabilistic uncertainty descriptions, depending on whether the uncertainty on the 
input is described using lower and upper bounds, or by a probability density function. The 
propagation approaches will be described first, followed by a brief overview of related topics, 
and finally a look at those techniques dealing with uncertainty amongst several models. 

2.1.1 Bounded Propagation Methods 

Bounded uncertainty techniques place lower and upper bounds on the model output (fig- 
ure 4). Gutierrez [11] describes several approaches. The first is the worst-case corners 
approach, in which the model is run at the lower and upper bound of each input parameter. 
Consider figure 5, where two parameters lead to four worst-case evaluations of the model. 
This is a relatively simple analysis with few parameters, though additional parameters 
quickly increase the computational burden. 

The first-order sensitivity approach in [11] uses sensitivity information between the 
output performances o and input parameters p to compute the performance bound based 
on small changes in p. 


do 

Once the sensitivities are computed, the performance uncertainty A o at different levels of 
A p can be determined. Since this approach relies on first order gradient information, its 
validity is only as good as that of the gradient. Often this limits the approach to small 
perturbations in p. 
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Figure 4: Uncertainty bounds provide a more complete picture than the nominal results 
alone; the left value nominally meets requirement, but with much greater uncertainty than 
the right value. 



P. 


Figure 5: Example of worst-case bounds, where a point is selected at the lower bound (L.B.) 
and upper bound (U.B.) of each parameter. 


Other techniques such as the Robust Control approach [11] and Linear Fractional Trans- 
formation [1] similarly produce a bounded uncertainty model. The key disadvantage of all 
of these approaches is that the bounds provide no information on the likelihood that a 
performance will be located at an extreme value. This makes it possible to reject a design 
that has an acceptably low but non-zero probability of failing to meet requirements. 

2.1.2 Probabilistic Propagation Methods 

A probabilistic uncertainty approach finds the probability density function (PDF) for model 
outputs, given the input PDFs. The uncertainty of an output can be described using 
standard probabilistic concepts such as mean and standard deviation. Instead of looking to 
see whether an upper bound exceeds performance requirement, the designer can report the 
probability that a design meets requirements. 

Monte Carlo sampling is perhaps the best-known propagation technique, in which ran- 
dom samples based on the input PDF are run through the model. Once enough output data 
points have been collected, their mean and moments can be taken to form the output PDF. 
Reference [10] provides a good overview of Monte Carlo sampling, along with other sam- 
pling techniques that try to reduce the large computational time. These include Stratified 
Monte Carlo sampling, Quasi-Monte Carlo sampling, and Latin Hypercube sampling. 

Melchers [18] describes the First- and Second-Order Reliability Methods (FORM and 
SORM), which determine the probability of failure of a system with normal random vari- 
able inputs. The limit state function, which defines the failure regions of the response, 
are approximated by either first or second order functions. Non-normal variables can be 
transformed to normal for approximate results. 

Linear Covariance Propagation (LCP) assumes that Higher Order Terms (H.O.T.) in 
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the design objective function are negligible: 


J{ a) = J 0 -f — + H.O.T. « Jo + ^ (2) 

where a is a vector of uncertain parameters. Under the assumption that a is a normally 
distributed vector random variable, a ~ iV(/x a ,E aQ ), the variance of the objective can be 
computed from the sensitivity, 


E[J) « J(m«) (3) 

. o-i dj „ aj r 

E[J ] « ^S QQ — (4) 

When the parameters axe independent, the contribution of each parameter to the vari- 
ance of the output can be computed: 


E\J 2 } ~ V ^ cr 2 
*=i OOLi 


da.i 


(5) 


where n p is the number of parameters. An example of linear covariance propagation in 
shown in Section 5. 

The change-of- variables routine involves the analytic PDF transformation [5] to prop- 
agate an input PDF through the model, producing an exact output PDF. This approach 
has been advanced in this research, and more information is provided in Section 4. 

Disadvantages of probabilistic uncertainty techniques include the computational time 
required for these approaches, especially for very large models. The other key problem is 
the accuracy of input distributions. While the Gaussian distributions is often chosen because 
of its simplicity, it may not be the appropriate distribution for a given input. Distributions 
are also difficult to determine empirically because of the large number of tests required for 
confidence in the distribution. Testing hardware in which only a limited number of copies 
are available may not be statistically sound. 


2.1.3 Related Fields 

There are several areas related to uncertainty analysis that also bear consideration. Robust 
Design [16], including the contributions of Taguchi, is a large field, concerned with methods 
of design and manufacturing to reduce defects. Stochastic Finite Elements [9] include 
stochasticity in the finite element formulation itself, although doesn’t appear to be as widely 
used as traditional finite element approaches coupled with uncertainty propagation tools. 

Design of Experiments (DOE) [20] is a field in statistics that deals with how to set-up, 
run, and analyze experiments involving random phenomenon. It is a well developed field, 
but rarely is included in any discussion of engineering uncertainty. Although originally 
developed for agricultural and biological testing, it has been brought into the realm of 
deterministic computational simulation to help examine stochastic behavior [10]. Computer- 
based sampling procedures such as Latin Hypercube Sampling are examples of modern 
DOE, where the goal is to gain the maximum information from a limited number of runs 
of a simulation. The high computational cost of running many of the integrated modelling 
simulations motivates further study of this field for uncertainty analysis. 
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2.1.4 Uncertainties between Models 


Most papers discussing uncertainty use one of the previous methods to consider the prop- 
agation of uncertainty through a single model. During a design lifecycle, many different 
models are used. New models are built with increasing complexity; data from testbeds are 
compared to models of the testbeds, and system components can be tested once available. 
As the models change, and as data becomes available, it becomes possible to map uncer- 
tainty across models and between similar structures. These require different tools than 
those used to propagate uncertainty through a single model. This branch of uncertainty 
analysis is much smaller, but several sources are available. 

Maghami [17] and Lim [15] have considered uncertainty between substructure models 
connected by component mode synthesis. They consider not only the propagation of uncer- 
tainty from one substructure to the next, but also the effect that the substructure assembly 
process has on the dynamics of the entire system. This is elaborated upon in the next 
section. 

Hasselman’s [12] database of aerospace structures is a compilation of knowledge gained 
from previous flexible structures for estimating the uncertainty of new models. The differ- 
ences between model and data from the historical structures are combined in covariance 
matrices. These covariances can then be used in a linear covariance propagation routine 
to estimate the similar differences in any new model of a similar structure. Bounds on the 
transfer functions of the new models axe produced. This approach was used by the SSL to 
examine model uncertainty for the Space Interferometry Mission [4]. The method does lack 
from the limited number of historical structures in the database. 

Campbell [6] describes uncertainty mapping between different configuration of a model. 
Using MACE as his example, he tests the structure on the ground and compares the results 
with the 1-g model. The differences are then mapped to the 0-g model in order to estimate 
the uncertainties of MACE as tested in the Space Shuttle. 

These three examples of using uncertainties between or across models are rarities in the 
literature. Most often only the single step of propagating uncertainty through a single model 
is shown. Any greater attempt at uncertainty management through the design process is 
performed on an ad-hoc basis. A common industrial approach, for which very few literature 
references have been found, is to use Model Uncertainty Factors (MUFs). These are little 
more than factors of safety about the nominal output. (See reference [2] for an example on 
determining this factor for dynamic modes.) 

2.2 Uncertainty Between Multiple Components 

Consideration was given on how to account for uncertainty between system components 
modeled separately. A quick description of methods to examine single component uncer- 
tainty is given first, followed by an overview of component mode synthesis to combine 
component models. Lastly these ideas are synthesized into a framework that could analysis 
the role of uncertainty between model. 
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2.2.1 Single Component Uncertainty 

Consider two components (figure 6). The dynamic response of each body, individually, is 
the time history y(t) given some forcing function F(t), which could be an impulse or entire 
time history. The response could be solved by integrating the dynamic equations of motion 


MiXi(t) + Ci±i(t) + KiXi(t) = /3iF(t ) 
yi(t) = Cy x x\ (t) 



Figure 6: Generic uncertain components with transfer functions from inputs F to outputs 

y 

for each body, i. The resulting dynamic response will be affected by any uncertainty in each 
model, either parametric in the M, C, or K matrices, or non-parametric in, for example, 
mismodeled elements or FEM fidelity errors. The input F(t) is often stochastic in nature; 
as a disturbance input it is usually modeled as unity white noise fed through disturbance 
shaping filters. Therefore, the standard measure of dynamic response is the root-mean 
squared value of the output. This can be found directly from the time histories by taking 
the expectation of the square of the output values. 

= y/E[y{t)*\ 

Instead of performing the integration, the RMS value can be found using the tools 
described in Gutierrez [11]. The easiest, perhaps, is the Lyapunov technique, for which the 
system must be described in state-space format 

x = Ax + Bd 
y — Cx 

where the input d is unity white noise, and the state vector x includes states from an 
appended white-noise shaping filter. Note that the i subscripts denoting the individual 
components have been dropped for convenience. The root mean square value of the dynamic 
response is found by solving for the state covariance matrix E 

A'Z + 'ZA 7 ' + BB t = 0 


and plugging it into 
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Gy = sfCYSJT 

A simple and direct method for determining the uncertainty of this value is to using 
random sampling (Monte Carlo). A sample is chosen randomly from a set of inputs with 
given probability densities, and after running a sufficient number of cases the statistics of 
the outputs translate into the PDF of the output system. The major drawback of this 
approach is the large number of test runs required to adequately cover the sample space. 

A more efficient method would be to use Latin Hypercube Sampling[10]. Each input 
probability density curve is subdivided into p regions of equal probability. If there are n 
uncertain parameters, then there exist p n bins. One approach, termed Stratified Monte 
Carlo sampling, would sample from each and every bin once. The actual location a point is 
taken from a bin could be a set point, such as the probabilistic center of that bin, or could 
be randomly placed. Stratified Monte Carlo sampling could still require a large number of 
runs, so to speed the process even further Latin Hypercube Sampling takes only one sample 
from each region of each input, and thus requires a total of only p samples. Figure 7 shows 
a graphical description of this. The results can be similarly analyzed to produce PDFs or 
CDFs of the uncertain output. 
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Figure 7: Illustration of Latin Hypercube Sampling 


2.2.2 Transfer Functions of Joined Components 

Either because different components are being build by different organization, or as a way 
to reduce the complexity of the global model, individual components are often modeled 
separately. In analyzing the system model, it is then often necessary to combine the com- 
ponent models in such a way that transfer functions G(s) = Y(s)/F(s) can be constructed 
where the input F is on one component, and the output y is on the other. 

There are several methods commonly used to assembly component models into a global 
system model. Blaurock [3] uses an acceleration feedback method, which inverts the inter- 
face inputs and outputs for one of the components, and treats the acceleration outputs from 
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Figure 8: Two structures joined 


the first model as force inputs to the second, and vice versa. Another popular method is 
Component Mode Synthesis [7]. This method will be described further. 

Component mode synthesis finds a basis of Ritz vectors to span the space of the interior 
and boundary (or interface) degrees of freedom for each component. The individual com- 
ponent states are transformed with these component modes matrices, and the transformed 
models combined by enforcing that boundary degrees of freedom are constrained to each 
other. 

The most popular variant of Component Model Synthesis is the Craig-Bampton method, 
first proposed by Craig and Bampton [13] with further elaboration in [14]. It obtains 
the transformation matrices by computing normal modes of the component interior points 
with the boundary degrees of freedom fixed. Constraint modes are found by setting each 
boundary degree of freedom to unity displacement or rotation, with no forces on the interior 
points. In the original papers the end goal of the analysis was to find modeshapes and 
frequencies of the unforced global system. Lim [15] considered substructure assembly from 
a controls sense, where inputs and outputs can be framed in a state-space manner. This 
was expanded upon in [17], where he also described uncertainties in the final synthesized 
models. These papers provide an excellent overview of the approach taken to find the 
assembled transfer functions. In addition the literature contains many other instances of 
component mode synthesis in use. 

A derivation of the input/output transfer function follows. It starts with the component 
models A and 5 in a generic form. 

Max a + Kaxa = PaF Mb'xbxb — PbF 

VA = CyAXa VB = CyBXB 

or, generally for any component i : 

Miii 4- KiX[ = PiF 

Vi ~ CyiXi 

Next, identify those degrees of freedom (DOFs) that will be attached to the other 
component. Separate the states into those at the boundary j and those in the interior i. 



then the main state equation matrices can be divided as shown: 

Mu Mij f X{ 1 ^ Ka Kij f 1 /3 p 

Mji 1 otj j K ji jfCjj ^ x j j 0 

Next write the physical coordinates in term of a Ritz approximation: 

Xi = ViPi 

where x are the physical coordinates, p are the generalized coordinates, and ^ are the 
pre-selected component modes. Methods of finding these Ritz vectors are described more 
fully in the literature on the subject, but generally use a combination of normal modes 
of each of the components, combined with a basis of constraint modes at the component 
interfaces. This differs slightly from the approach in [15], which includes as an input the 
boundary forces due to the next component, and finds the proper “control loop” such that 
the boundaries are constrained together. 

Once the redundant coordinates are removing using a second state transformation p = 
Sq , the end result is the state-space representation of the assembled system, 



The system can also be written in transfer function form. If the input is on component A , 
and the output on component B , then (3 a has a value of 1 at the degree of freedom of the 
input and zeros elsewhere, and C y s has a 1 at the degree of freedom of the output, and 
zeros elsewhere. (3b and C v a are all zero, in this case. The transfer function is: 

G AB (s) = ^^=[ 0 C yB } ■* b .S.[s 2 M s + sC s + K s ]- 1 -S t -^- ^ 

The inverse transfer function, with input on B and output on A is: 

g ba(s) = ^^-=[0 C yA ] * A -S [s 2 M s + sC s + K s ]- 1 -S T -* T B - Pg 

And since C v a = /3^, and C y s = 0$, these transfer functions are also the transposes of 
each other. The transfer functions should be the same, although some variations due to 
removed normal modes are likely. 
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(a) Joining of generic components 


(b) Interface of MACE 


Figure 9: Component interfaces 


2.2.3 Interface Uncertainties 

Next the modeling uncertainties at an interface of two generic structures are considered. On 
a more specific level, one could imagine the main boom of MACE being modeled separate 
from the gimbal and flexible appendage. The model would be assembled (using the Compo- 
nent Mode Synthesis approach proposed, for example) with the interface at the node- gimbal 
attachment degrees of freedom. 

Along with all of the typical uncertainties in the models themselves (parametric, model 
error, etc.) uncertainties in the global system could also result from the assembly process 
and the interface itself. 

The first source of interface uncertainties could be from the component mode synthesis 
method. One of the advantages of the method is that it reduces the number of states of 
the component models by keeping only a small subset of the normal modes. The final 
assembled model can then be constructed with only the critical states of each component 
retained. Although this reduction capability is an asset in working with large models, it 
has the potential to allow important dynamics to be lost in the reduction. Maghami and 
Lim [17] suggest that component modes up to twice the desired system bandwidth should 
be kept, but admit that even this will not guarantee accuracy in all cases. The uncertainty 
resulting from these eliminated modes must be considered. The easiest solution would be 
to run convergence tests, with different numbers of kept modes. Once the results begin to 
diverge, the acceptable limit of mode reduction has been reached. 

A similar procedure could also be used to determine the effect of the choosing the 
wrong modes. The Craig-Bampton method uses normal modes with fixed boundaries and 
constraint modes, but there are other options. Normal modes can be found with fixed, free, 
or intermediate boundaries. Free body modes must be accounted for in the case of free 
body motion. Reference [7] describes some of these other options. 

Another source of uncertainty could be model errors in the interface itself. Although 
the illustration in figure 9(a) shows the interface as merely a single point, physical inter- 


faces could have static and dynamic behaviors of their own. For example an optics on a 
telescope may be separated from the noisy ground by an isolator assembly with its own 
dynamics. Thus, depending on the complexity of the interface, some or all of the standard 
model uncertainties might need to be included, including parametric and non-par ametric 
uncertainties. 


interface 



Figure 10: Component interface 

A simple model of an interface would include only stiffness terms, and no dynamics. 
The standard illustration could be redrawn to show this (figure 10). The interface would 
become a part of the assembled system using 

A K r A 

M = 0 K = Ki 

MrB MrB 

Given Interior and Boundary states [ x\ xf ] T for each component, the components 
can be assembled using the transformation 

J 0 0 " 

0/0 
0/0 
0/0 
0 0 I 

The interface stiffness uncertainty could then be described as K\ — K 4- A K for use in 
uncertainty calculations. 

2.2.4 System Performance Uncertainty 

Consider two alternate models for the assembled system. The first is a networked model in 
which all of the component models are fully retained as distinct entities. 

This is akin to the integrated modeling approach, where complete, distinct models are 
combined together to form the total system response. Since the full system model for each 
component remains, all of the parametric uncertainties of the individual models can be 
retained and used to find the total system performance uncertainty. Both of the approaches 
described above are utilized in this case. The models axe combined using Component Mode 
Synthesis. In addition to the uncertainties in the models themselves, interface uncertainties 
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can also be included. The combined model in state-space form gives the system performance 
using the Lyapunov disturbance to performance algorithm. 

Given a collection of probabilistic uncertain parameters, the Latin Hypercube Sampling 
technique described initially can be used to examine the uncertainty space much faster 
than traditional Monte Carlo sampling. More parameters can be investigated than would 
be possible with the Change-of- Variables routine. One could also consider using First or 
Second Order Reliability Methods [18] (FORM/SORM). 

The other approach to take is to build an aggregate model of the system which does not 
retain the individual physical parameters of the component models. This “tree” method 
of modeling, would usually result from models being built by different divisions or or- 
ganizations. The separate models are then independently reduced from the full physical 
coordinates (thousands to hundred of thousands of degrees of freedom, perhaps) to a re- 
duced set of non-physical coordinates. This reduction could use component mode synthesis, 
although other techniques such as Guyan reduction are also possible. In any effect, the re- 
duced models are assembled into a global systems model elsewhere (up the ’’branches” of 
the tree-like organization, say) without reusing any of the lower-level component models. It 
is common in this case that reduced mass and stiffness matrices ( M r ^ and K r B , using the 
notation above) are passed to the assembly process. 
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Figure 11: Tree model 

In order to compute the total system uncertainty, the individual uncertainties in the 
components must still be passed along with the mass and stiffness models. It is therefore 
required that the parties responsible for the component models also prepare an uncertainty 
model. 

One of the challenges in this case is that the physical coordinates are no longer available 
for the aggregate model. Therefore, the physical uncertainties used previously cannot be 
passed along. The challenge is to frame the uncertainties in the transformed coordinates, 
and pass these along as A M t a, A K t a terms. As shown in [11], it is relatively easy to use 
transformed, modal uncertainties such as modal mass, frequency, and damping, as uncertain 
parameters. In terms of component mode synthesis, [15] and [17] both choose to use only 
normal modes for the component mode basis. In this case, the reduced systems are in the 
familiar form 
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where m r i is the modal mass, Q is the damping, and t Qi is the natural frequency for com- 
ponent i. Using the transfer function from the component input to either the component 
output or the interface, the largest modal contributors to the cumulative RMS should be 
identified by the component organization, and singled out for inclusion in the component 
model. 

In order to consider probabilistic uncertainty, it would be necessary also for the compo- 
nent organization to include information on the probability density of the uncertain modal 
parameters. This could be considered in two ways. Given that the component organization 
has access to all of the uncertain physical parameters, a probabilistic propagation routine 
could be used to directly compute distributions for critical modal parameters. These dis- 
tributions could be passed along with the models. If there are a large number of uncertain 
parameters passed along however, it could be reasonable to invoke the central limit theorem, 
and assume that the modal parameters would tend toward a normal (Gaussian) distribu- 
tion. Only the mean value for each modal parameter, along with a variance estimate would 
have to passed along in this case. 

The model handed up the tree to the system assembly division would be of the form 


M ri = M ri + A M ri K ri = Kri + A Kri 


for component i, where the M ri and K r i are mean values, and the distribution is described 
in the A terms. As an example, if the frequency were normally distributed about a mean 
value contained in K r u with standard deviation the uncertainty would look like: 


A K ri = 


0 

N(0 ,<7^) 


0 


The aggregate model can now be constructed without use of the lower level component 
models, still by using the assembly procedure described earlier for component mode syn- 
thesis. The performance can be computed using the same Lyapunov techniques described 
earlier. The same approach as before can be taken to compute the uncertain performance, 
namely by use of Latin Hypercube sampling over the uncertain parameter space. The key 
difference in this tree method is that the uncertain parameters are now in a transformed 
coordinate space which is much smaller than the original physical coordinate space. 


2.2.5 Conclusion 

Two methods of computing the performance uncertainty of a system comprised of two or 
more uncertain components was presented. The first method, the networked method, uses 
the entire component model and has at its discretion all of the uncertain physical inputs. 
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This is both an advantage and added difficulty. The advantage is that the global sys- 
tem modeler has a great amount of insight into the actual physical causes of system-level 
uncertainty. The design space has much greater fidelity since the effect of individual uncer- 
tainties cam be identified. Resources cam be brought to bear immediately on the physical 
sub-component responsible for the system uncertainty. Also, there is less approximation in 
this analysis technique, as compared with the tree method which must first approximate 
the uncertainty in terms of modal parameters. 

These advantages of course are tempered by the great size of the full component models 
being used. Running the analysis will take more computer power. Also, there is the very 
real possibility of being overwhelmed by the sheer number of possible uncertain physical 
parameters. The critical effect of one component on another very well could be reduced to 
several important modal frequencies and a full physical system may not be necessary. 

Model assembly using the tree method starts with the models reduced already. Uncertain 
modal parameters must be decided upon and included with the reduced models. The 
advantages here are that much smaller models are used in the uncertainty analysis, so 
that either less computer power is needed, or more probabilistic samples can be run. The 
system modeler also has a much simpler task of assigning uncertainty, since he/she is not 
responsible for the details of the component, but rather of only the critical effects. The 
disadvantage is of course loss of insight into causes of performance uncertainty. The system 
modeler can only isolate the component which creates the most uncertainty. Identification 
of the actual input causing the component to be uncertain can only be done by sending 
word back down the tree to the component group, which must investigate the problem and 
send back updated models. This process can take much longer than the networked method. 

2.3 Role of MACE in Uncertainty Analysis 

The MACE testbed is used extensively this research. It provides an excellent opportunity 
to examine uncertainty on an accelerated timeline, since the hardware is already designed, 
developed, and tested. Data are available in both 1-g and 0-g environments, across many 
years and with different sets of hardware. The presence of this testbed is a great advantage 
in uncertainty research, since any validation of uncertainty routines requires knowing the 
actual behavior of a modeled system. This access to working, tested hardware allows MACE 
to be “designed” anew, with immediate comparisons between those new designs and actual 
behavior. The design process can be recreated without the years of manufacturing normally 
required. In order to prevent knowledge already gained from MACE from influencing the 
uncertainty models of the new design, actual flight data will be kept separate until such 
time as would be appropriate in the design process. 

The design process has been started anew on MACE. A new model has been constructed 
using Femap® (figure 12). It is based solely on original drawings of the MACE article, 
and uses handbook values for all materials. This model is akin to a preliminary “stick” 
model, often used to perform the very first dynamical analysis on a structure. There 
also exist MACE models of greater accuracy, whose parameters have been updated with 
data. While these models may be more accurate than this preliminary model, they also 
rely on data from hardware which may not exist in a real design process. These models 
together represent several stages of the the design process and are a valuable tool in testing 


17 



uncertainty approaches and in tracking uncertainty through a design. The updated MACE 
model is used in Section 5 to illustrate uncertainty analysis based on hardware testing. 



Figure 12: New finite element model of MACE 


3 Uncertainty Source Investigation 

The first step in placing uncertainty bars on the performance output is to identify those 
aspects of the model inputs that are most uncertain. Generally these may include both 
parametric uncertainties, representing all parameters explicitly contained in the model (e.g. 
Young’s modulus or rotary inertias) and non-parametric uncertainties, defined as all other 
uncertainties that are not explicitly in the model (non-linear behavior, incorrect boundary 
conditions, etc.). 

A list of many typical sources of model uncertainty in high performance optical systems 
was produced (table 1). It was made from a combination of literature on the subject and 
SSL experience dealing with uncertainty. The intention was to produce a comprehensive 
list which could be compared with a new design model; those items on the list which are 
present in the system can be given further consideration. 

While this list identifies many possible sources of uncertainty, it is still necessary to 
identify those particular sources most likely to affect a particular design. A goal of un- 
certainty research should be to create a quantitative tool that can parse through such a 
comprehensive list and highlight sources for a model. Currently, engineering judgement is 
still necessary to decide which sources to investigate in greater detail. 

Using MACE as an example, several uncertainty sources were chosen based on their 
expected influence on the model (all are expected to contribute to the final performance) 
and by the fact that all are relatively uncertain. These parameters are given in table 2, and 
consist both of parametric uncertainties such as modulus, and non-parametric uncertainties 
such as non-linear behavior in the strut-node connection. 

One of the greatest difficulties with uncertainty modeling is that even when sources of 
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Table 1 : Uncertainty Sources 


Material 

Modulus 

Density 

Point mass 

Damping coefficient 

Composite Ply Orientation 

Composite matrix chemistry 

Model Error 

Mismodeled loading 
Modeling simplication 
Rigid links 
Rigid substructures 
Beam approx, for truss 
Element formulation 
Boundary conditions 
Cross-inertias (bending, rotary) 
Sensor misalignment 
Rotation, translation 
Interface compliance 
Preload 
(de)stiffening 
Static deformation 
Eccentric loading 
Joint rigidity 
Gross incompetence 


Gravity 

Joint Pre-load (locking) 
Body forces 
P re-stiffening 
Sag (coupling) 
Gravity-sensor coupling 
Suspension dynamics 
Pendulum modes 
B/C compliance 
B/C damping 
B/C energy leakage 
B/C noise sources 

Degradation 

Fastener hole wear 
Change in modulus 
Cracks/crack density 
Bearing wear 

Manufacturing/Testing 

Geometric tolerance 
Material residual stress 
Fastener Torque 
Preload 


Environmental 

Temperature varying E 
Temperature varying v 
Temperature varying ^31 
Preload (thermal strain) 
Moisture (mass, shape change) 
Outgassing 
Air dissipation 
Density variation of air 
Speed variation of fight in air 
Acoustic noise 

Non-Linear 
Modulus non-linearity 
Bearing stiction 
Bearing rattle 
Bearing imbalance 
Bearing shape irregularity 
Joint stick/slip 
Loose fastener 

Policy 

Cost 

Miscellaneous 

Optical effect of air 


T able 2: Identified Sources of Uncertainty for MAC E 
Parameter NASTRAN card 

Young’s modulus MAT1 

Rotational stiffness CELAS2 

1 12 cross-bending inertia PBAR 

Non-linear interface 
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uncertainty are identified, very often there are no accurate models of how the source varies. 
This problem is particularly acute with probabilistic descriptions of uncertainty, as there 
are few references that give PDFs for common model inputs such as modulus. This lack 
of accurate uncertainty inputs casts a shadow on all of the analysis results, since without 
physically meaningful inputs the outputs are of dubious value. 

Therefore, one of the goals of this work is to attempt to experimentally measure the 
distribution of prominent uncertainty sources. Two methods were used: direct testing of a 
MACE strut and model updating of selected NASTRAN input parameters. 

3.1 Direct Measurement 

If possible uncertainties are identified and the hardware exists, a model of the variation of 
the parameter or effect of the non-parametric error can be produced by empirical testing. 
This has been done with a node-strut-node component of MACE. The first source of un- 
certainty under investigation is the Young’s modulus E of the polycarbonate making up 
the length of the strut. Modulus was used in the model updating routines performed with 
the MACE model and allowed very good fits of the model to the data. Despite this, it was 
not clear whether the updating routine pushed “too hard.” and thus varied the value for 
this parameter outside of the physically realizable range. The other source of uncertainty 
investigated is possible non-linear behavior at the strut-node interface. Depending on the 
tightness of the the steel collar that holds the aluminum node to the polycarbonate strut 
(figure 13), it is expected that slip between the strut and node could manifest into non-linear 
behavior. 



Figure 13: One of the MACE nodes, with two struts and collars connected 

The first tests performed were axial compression tests using an Instron testing machine 
(figure 14). Since polycarbonate is an isotropic material, no difference in modulus was 
expected between bending tests and compression tests. The Instron machine applied in- 
creasing load at a rate of two pounds per minute. The load varied between -10 and -110 
pounds. Load and displacement were measured by the testing machine, and strain was 
measured by a strain gage attached on the center of the strut. Example stress-strain and 
force-displacement plots are shown in figure 15. The stress was computed using the cross 
sectional area of the strut 


( 7 = 


P measured 


7 rr 


2 — r 2\ 


( 6 ) 
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Figure 14: MACE strut compression test 
where r D is the outer radius of the strut, and r* the inner radius. 



Strain jjt Displacement pr] 


(a) Stress-Strain test (b) Force-Displacement test 

Figure 15: MACE strut compression test results 

Upon inspection of the plots, it is obvious that the displacement data is much more 
noisy than the strain gage data, but trends can be seen in both. There is hysteresis which 
indicates non-linear behavior; this is not, however, considered to be the node-strut non- 
linearity expected. The interface slip should result when the load goes through zero: since 
only compression is being applied, that never occurs. In this case, the hysteresis in the force- 
displacement plot results from the strut being pushed axially into the node in compression. 
Since the strain gage measures strain just on the strut, the amount of hysteresis is much 
smaller than the displacement sensor which measures absolute displacement of the entire 
assembly. Because the load never goes through zero however, this test is only useful for 
measuring modulus. 

By measuring the slopes of the resulting stress-strain curve, the averaged value for the 
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modulus was 117,000 psi. The data sheet for the particular material used (Zelux polycar- 
bonate by the Westlake Plastics Company) lists a flexural modulus of 340,000 psi, nearly 
three times the value measured. This gross difference was extremely suspect, and was not 
taken as physically meaningful. Especially since the primary loads on MACE would be 
bending, rather than compression loads, it was decided that the test was incorrect for the 
entire component. 

Following the failure of the compression tests, a bending test more akin to the expected 
loading on MACE was set up (figure 16). The strut w T as clamped to a laboratory bench. 
A L-beam was placed between the table and the fixed node in order to provide greater 
boundary condition rigidity. A point load was applied at the opposite node in order to 
produce bending. The load was applied by a stinger attached to a shaker. A loadcell was 
placed in the load path to measure resulting force. 



Figure 16: Strut bending test setup 

The load was applied at 0.1 Hz and in such a manner that the displacement of the end 
node would go through zero. Bending strain was measured using strain gages at the center 
of the strut, and the stress is computed using Bernoulli-Euler beam theory with a point 
load. 


V 

^x(^) = j [~ ft stinger T ( L — X )] Ploadcell (7) 

*zz 

r Q is the strut outer radius, I zz is the area moment of inertia, L is the length, Ploadcell is 
the force from the load cell, and S st i n ger is a corrective term that takes the stiffness of the 
stinger into account. 

The struts used are those of the MACE Engineering Model (EM), which was built for 
hardware testing before construction of the flight hardware. Three EM struts were available 
for testing, allowing difference between otherwise identical components to be examined. In 
order to determine the effect that collar tightness has on any non-linear behavior, the 
collars were turned to three levels of tightness, qualitatively defined as 'loose' (threads 
barely engaged), 'normal' (collar snug against the node), and Tight 7 (as tight as possible by 
hand). Data was taken over multiple days, over a two month period. 

An example stress strain curve is shown in figure 17. The most immediate observation is 
that the system is highly linear for all levels of collar tightness. Further, the collar tightness 
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seems to have no appreciable effect on the results. Even as the stress goes through zero (as 
the strut changes from bending downward to upward or vice-versa) there is no non-linearity 
detected. This is an important result, since it eliminates one of the potential sources of 
uncertainty. Now that non-linear behavior in the strut-node interface is shown not to be a 
factor, resources can be brought to bear examining uncertainties elsewhere. 



Figure 17: Stress-strain results from bending tests. EM strut #1, at three levels of collar 
tightness. 

Modulus data is also taken from the slopes of the stress-strain curves. All data collected 
is shown in figure 18. Each strut is identified by a color. The collar tightness differences 
are again seen to be negligible. The modulus values are much closer to the published 
value, and taken as physically significant. The experimental moduli are all higher than the 
published value, between 375,000 psi to 403,000 psi. Figure 18 more significantly shows 
the differences between the three struts. The variation within a single strut is never more 
than 7,000 psi, while between struts can be upwards of 22,000 psi. This suggests that the 
subtle differences between otherwise identical hardware dominate the uncertainty. These 
differences may result from manufacturing tolerances or from the fit of the strut against the 
node. Likewise, a possible source of the difference could be sensor differences in the strain 
gages, and not a physical variation at all. This illustrates the difficulty in isolating the true 
source of component uncertainty. 

The statistical mean and standard deviations of the individual strut data, along with 
the aggregate data, is given in table 3. The statistical significance of this can be questioned 
because of the very limited number of struts available. Whereas material coupon testing 
could use hundreds of samples to give confidence in the statistics of material modulus, tests 
of large components often have to make due with a limited number of hardware pieces. 
This restriction means that it would be difficult to produce detailed probability density 
functions with such samples. A much larger uncertainty effort, aimed at creating a true 
statistical database of material properties and other parametric model inputs would be 
highly advantageous to the field. 
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Figure 18: Averaged strut moduli for the three collar tightness levels. Colors represent 
strut #1 (blue), #2 (green), and #3 (red). 


Table 3: Strut bending data statistics 



Mean 

Std. Deviation a 

Strut #1 

377,000 psi 

1700 psi 

Strut #2 

386,000 psi 

2500 psi 

Strut #3 

399,000 psi 

1600 psi 

Aggregate data 

387,000 psi 

9300 psi 


Given the limited sample size with only the mean and standard deviations of the data, a 
Gaussian PDF is an acceptable distribution. Figure 19 shows the probability mass functions 
(PMF) of the discrete data points, overplotted by Gaussian PDFs. 

The results of the component uncertainty tests show how difficult it can be to isolate 
and produce accurate distributions of identified uncertainties. However once uncertainties 
such as modulus or non-linearities are suspected, testing proves invaluable. The lack of 
non-linear behavior in the interface means that no more effort must be wasted in examining 
it. Even with the questions that remain on the statistics of the modulus data, the results 
show that the effective modulus is higher than the published values. A Gaussian PDF or 
3<t bounds on the mean values still provides uncertainty information that would otherwise 
be complete guesswork. 

Lessons learned include the need to create tests that properly capture the uncertainty 
under investigation. The MACE struts are loaded primarily in bending, so a bending test 
should have been performed first. For probabilistic analysis it is also necessary to have a 
large enough sample size to ensure statistical significance. For tests of larger components 
for which only a limited number of pieces are built, probabilistic analysis may be limited to 
Gaussian PDFs. 
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Figure 19: Gaussian PDFs from the bending test results with the individual struts in colors, 
and the aggregate data in black. 

3.2 Physical Parameter Fitting 

The second approach used to build uncertainty models is to use frequency response data to 
update model parameters via a tuning methodology. Using data from multiple test runs, a 
collection of parameter fits can be used to create a statistical model of those parameters. 

The parameter tuning methodology in general is an iterative process, with a “coarse” 
stage in which large model changes are made to get the gross features of the model to 
match the measured response of the system. Then a “fine” stage is used to fine-tune the 
model. First, the tuning metric is chosen. Modal frequencies and Frequency Response 
Functions (FRFs) are good candidates for the coarse and fine stages, respectively, of the 
tuning process. Next, tuning parameters are chosen using engineering insight, and the 
functional dependence of the tuning cost on the parameters established. This stage is 
itself iterative, to find parameters that have the necessary effect on model response (modal 
frequencies and residues for example). Next the parameters are tuned to create a “best- fit” 
model, either manually or via nonlinear search methods. If the best-fit model is sufficiently 
accurate, the process is stopped. If not, the parameter set is not complete and additional 
tuning parameters are identified. 

The parameter model is generated by creating multiple tuned parameter sets, once the 
complete set has been established. 

The parameter model for the MACE MBP was created from physical parameter fits 
to three different data sets, taken on three different days after complete disassembly and 
reassembly of the MACE hardware. The parameters were chosen as those with the largest 
influence (highest sensitivity) in the FRF channels. The mean and standard deviation 
of the fits were computed using MATLAB standard functions, and independent Gaussian 
models were fitted to each parameter. The parameters are tabulated in table 4. The 
identified parameters and 3-sigma bounds are shown in figure 20. The parameters are: 
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Physical parameter fits and 3 sigma bounds around mean 


Active skill piezo - 
App 112 
PrAppX 
Pr App Z " 
Sc App X ' 
Sc App Z * 


Data Set 1 
Data Set 2 
Data Set 3 


Figure 20: Physical parameter fits (plotted as deviations from mean of identified value) and 
3-sigma bounds of Gaussian model. 


Lexan end modulus, that models the node-to-Lexan strut interface; Lexan center modulus 
that models the struts themselves; the active strut piezo modulus; cross-inertia in the 
appendages that cross-couples bending motions; and Primary and Secondary Appendage 
base rotational stiffnesses that model compliance in the attachment between the appendages 
and the gimbals. Note in particular that the Primary Appendage base X stiffness has a 
very wide distribution. The variation in the fits was driven by the need to match a highly 
observable mode that changed by a few Hertz in frequency from test to test, and that was 
traced to a machining tolerance that left the appendage attachment screw slightly out of 
flush with the face of the gimbal. The result was a preload that was highly dependent on 
the tightness of the screw. 


Table 4: Physical model update parameters 


Name 
Lexan ends 
Lexan center 
Active strut piezo 
App. I 12 
Pr. App. X 
Pr. App. Z 
Sc. App. X 
Sc. App. Z 


Description 

Lexan/node attach point modulus 
Lexan strut barrel modulus 
Weighted piezo strut modulus 
App. bending cross-inertia 
Primary App base X rotational stiffness 
Primary App base Z rotational stiffness 
Secondary App base X rotational stiffness 
Secondary App base Z rotational stiffness 


Mean 

2.51 GPa 
2.77 GPa 

3.52 GPa 
53.4e-12 m 4 

444.3 Nm/rad 
762.0 Nm/rad 
731.5 Nm/rad 
832.2 Nm/rad 


Std. Dev. 
0.25 GPa 
0.11 GPa 
0.48 GPa 
4.90e-12 m 4 
84.1 Nm/rad 

12.7 Nm/rad 
59.9 Nm/rad 

33.7 Nm/rad 
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4 Change-of- Variables Propagation Routine 

The Change of Variables solution uses the Probability Density Function (PDF) transfor- 
mation [5] 

/ £ ( ct ) = fp(p) (8) 

to describe the PDF of the output u to inputs p, given the PDF of the inputs, and assuming 
a relationship a = g{p). Generally this requires the functional relationship to be inverted 
p = h(a), so that 

M<r) = \^\fp[h(cr)) (9) 

For systems analysis, a is often a performance found using frequency domain methods, 
with the inputs p buried in the process; in such situations is can be difficult to explicitly 
form the inverted relationship p = h(a) and the sensitivity dp/da. However, when the 
function is known implicitly, the form can still be used numerically to sample the PDF of 
a . The approach is to rewrite the differential relationship as 

• h{a) = fpip) (10) 

This is an algebraic relationship that can be solved for /s: 



/s( CT ) = 


fpjp) 

da 

"Sr 


( 11 ) 


which only requires the computation of the gradient of a with respect to p. In the case of 
multiple input-multiple output systems, this becomes the determinant of the Jacobian J. 


OPl ’ ' * &Pn 

J = det : (12) 

dcr n d(T n 

- dpi dp n - 


Note that an equal number of inputs and outputs are required so that the Jacobian remains 
square. If there are smaller number of output, “dummy” outputs must be included at this 
stage, and can be integrated out of the results later. 

This approach to an output distribution is often described in textbooks on the subject, 
but is not found in research literature. Examples most often include simple single input- 
single output functions that can be solved analytically, with pencil and paper. This research 
has taken this procedure and applied it to much larger problems that cannot be solved as 
easily. An important step is to sample the functions, instead of attempting to perform the 
inversion. Given a range of values about the inputs p, the output cr, the input PDF fp(p) 
and the determinant of the Jacobian J can be computed at each value, equation 11 is used 
to compute the output PDF at each point. Finally, the output PDF values are mapped from 
the input p-space to the output cr-space, where they must be binned into the appropriate 
output grid space. The exact procedure is described next, using inputs x and outputs y for 
a relationship y = g{x). 
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1. Grid input x 

2. Pick x Q 

3. Compute fx( x o) 

4. Compute y Q = g(x Q ) 

5. Compute |^| 

6. Compute fy(y 0 ) - fx(x a ) • 

7. Repeat steps 2-6 through the x-grid 

8. Grid output y 

9. Find all yi in each Ay bin 

10. Sum together all fyiVi) that correspond to those yi in the bin for the final output 
distribution 

4.1 Monotonic Example 

This procedure will be demonstrated first with two simple examples. A full parameter-input 
to performance-output dynamic analysis will be shown in Section 5. First consider a system 
of equations described below, in which the inputs and outputs are related monotonically, 
i.e. through the sample space the gradients are never zero. 


yi = xi - x 2 (13) 

y 2 = 4a:i + xix 2 


This procedure makes no assumption on the distributions of the input parameters. For 
these examples, the Gaussian distribution is used to model the input PDFs. 


fxAxi) = 



( -foi ~Mx,) 2 \ 

l ) 


(14) 


If the inputs are independent, their joint PDF can be formed by multiplication. The result- 
ing distribution is plotted in figure 21. 


fx uX 2 ( x U x 2) = fx i(*l) • fx a (x 2 ) (15) 

The gradients can be found using the sensitivity functions of the DOCS toolsuite. For this 
simple example, the determinant of the gradient is quickly calculated analytically. 


Jana, — X 1 H - x 2 d~4 

And finally the output distributions in the x-space are computed via 


(16) 
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Figure 21: Joint PDF for two Gaussian inputs 


fY*y*(yi,V 2 ) 


f X 1 ,X2 (*^1 i *^ 2 ) 

PI 


(17) 


This output PDF is shown in figure 22(a). Once the PDF values are mapped from the 
input-space to the output y-space, the PDF looks like figure 22(b). This surface is the 
joint PDF for both outputs. In order to solve for the individual PDFs f\\ and fy 2 (if, for 
example, 1/2 was a •‘dummy” output), the undesired output can be integrated out. 


fvdyi) = j_ 


fY 1 y 2 (y\-y2) dy 2 


(18) 


This procedure can be done using numerically integration. The resulting individual PDFs 
and associated cumulative distribution functions (CDF) are shown in figure 23. Both CDFs 
asymptote to 1.0, guaranteeing that the output distributions are legitimate PDFs. 



(a) Output PDF, whose points are mapped to (b) Output PDF, whose points are mapped to 
the input grid the output grid 


Figure 22: Output PDF mapping 



Figure 23: Output probability density functions (PDF) and cumulative distribution func- 
tions (CDF) for the monotonic example. 
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4.2 Non-Monotonic Example 

The second example is of a similarly simples set of equations, but in this case the equations 
are non-monotonic. This results in gradients that go to zero. This will cause difficulties in 
equation 11, with J in the denominator. 

The example system of equations is 



(a) Function yi(x\,X2) (b) Function 1/2 (^1,^2) 

Figure 24: Example system of equations. 

plotted in figure24, with analytical J as below. A plot of J is provided in figure 25. 

J = 2x\ - 4x2 — 2x2 (20) 

The difficulty is that the surface J (multi-dimensional if there are more than 2 inputs/outputs) 
goes through zero. As the steps are performed to compute the output PDF, a line of singu- 
larities appears. This can be seen as the PDF is plotted against the x points in figure 26(a), 
and again against the y points in figure 26(b). At this line the output PDF “wraps” around 
itself, creating two surfaces which must be summed together for the final output. It is 
also necessary to create a much finer mesh around the singularities, in order to capture the 
complete PDF; a coarse mesh will produce much higher values than actually exist. The 
Change-of- Variables code built allows for variable meshing, and provides the binning capa- 
bility for non-monotonic systems such as this. The final resulting PDF for this example is 
shown in figure 27. An example using the MACE hardware is shown next. 
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*1 % 


Figure 25: Determinant of the Jacobian. Note that in the non-monotonic example, a line 
of the surface passes through zero. 



(a) Output PDF gridded in x-space (b) Output PDF gridded in y-space, before 

binning. Note that the surface from plot (a) 
is wrapped about itself in this grid system. 

Figure 26: Output PDF samples 
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Figure 27: Resulting output PDF from CoVJ2D.m. 
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5 Application of Uncertainty Analysis in the Design Process 

All of the tools outlined in this report: legacy propagation tools, measured uncertainty 
parameter distributions, and the change-of-variables propagation routine, were used to il- 
lustrate an uncertainty analysis methodology with the MACE hardware. An overview of 
the framework developed in shown in figure 28. 



Critical 

parameter 

selection 


Critical 

parameter 

propagation 


Figure 28: Uncertainty analysis framework. Data gained from component testing or hard- 
ware FRFs are used to generate a statistical model of dynamic model parameters. Resulting 
parameter sets are statistically analyzed to generate a parameter model. Linear Covari- 
ance Propagation is used to find the largest contributors to the performance variance, and 
Change of Variables analysis is used to determine the exact output PDF and CDF. Results 
are checked using Monte Carlo simulation. 

The uncertainty analysis procedure used in this example is: 

1 . Conduct multiple experiments under varying test conditions (time, temperature, dis- 
assembly/assembly) to create multiple datasets of nominally identical hardware. 

2. Perform a physical parameter identification on each dataset, in which parameters of 
the physical model are tuned to best fit the data. Create a statistical model of the 
identified parameters (by choosing the form of the distribution, and coefficients). 

3. Perform an uncertainty contribution analysis using Linear Covariance Propagation 
(LCP) to identify the parameters that create the largest variance of the system design 
objective function. 

4. Perform a Change-of- Variables (CoV) analysis on the parameters with the largest 
contributions to determine the exact design metric distribution. 
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5.1 Test System 

The design approach is evaluated on the MACE Multibody Platform (MBP) in zero-gravity. 
Parameter studies and model updating were performed as part of the MACE Reflight, and 
provide the parameter database. This is a useful sample problem because it exhibits many 
of the features important in high performance controlled structures: 

• coupled structural and control system dynamics, via sensor and actuator dynamics, 
and computer time delays 

• closed loop control of the servoed gimbals and reaction wheel actuators 

• high bandwidth flexible mode observability in the performance 

• realistic scaling issues in the form of a wide frequency band (less than 1 Hz to over 
100 Hz) and modal observability (Hankel singular values ranging from 495.8 to 8.9e- 
10 ) 

The problem demonstrates propagation of physical uncertainties through a balanced 
state reduction and transformation to real-modal form. Additionally, the system will even- 
tually support closed loop high bandwidth control studies. 

The parameter model is that described in Section 3.2. The physical model parameters 
chosen by engineering insight, and updated using model tuning, are modulus of the lexan 
center and ends, modulus of the active strut piezo, appendage bending cross-interia, and 
rotational stiffnesses of the appendages. Table 4 gives the mean and standard deviations 
for all of the parameters. 

5.2 Objective Function 

The design objective for the MACE mission is to hold the 2-axis inertial pointing angle 
of the Primary Appendage (PrApp) under 2-axis white noise commanded slewing of the 
Secondary Appendage (ScApp). This is posed as the Root Sum Squared (RSS) PrApp 
pointing error as measured by the Primary Rate Gyro X (PRX) and Primary Rate Gyro Z 
(PRZ) sensors, which is a scalar cost: 

J = yjE[PRX 2 } + E[PRZ 2 } (21) 

The scalar RMS of PRX and PRZ will also be examined to understand the relative 
importance in the cost. 

The units of the design model are nanometers to degrees. For the purposes of evaluating 
probability of success, the RSS requirement is set to be 5% larger than the predicted RSS 
(in other words, that the system achieves the requirement with a margin of 5%). 

A unit torque intensity is assumed, which is within the capability of the actuators. In 
order to avoid exciting the low-frequency dynamics (due to the tethers which “stationkeep” 
the MBP inside the Middeck), the disturbance actuators are shaped with a filter that 
removes low (below 5 Hz) and high (above 200 Hz) frequency disturbances (figure 29). This 
has the effect of making the flexible mode dynamics more important in the pointing RMS. 
Without the filter, low-frequency components in the disturbance cause excessive pointing 
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errors due to rigid body motion, which wash out any flexible dynamics which are functions 
of the structural parameter uncertainties. 


1 Mm RMS Disturbance shaping filter 



Phase 



Figure 29: Disturbance shaping filter used to de-emphasize low frequency modes. 


5.3 LCP Results 

The first step is to evaluate the combined effect of all uncertainties in the quick but approx- 
imate LCP analysis. 

The CDF and PDF of the RSS pointing error, computed using LCP, is plotted in 
figure 30 (top and bottom respectively). The horizontal red line in the top plot shows the 
probability of meeting the requirement of 5.99 deg (vertical red line) which is 79% in this 
case. The vertical red line shows the mean output RSS pointing error of 5.70 degrees, with 
variance of 0.35 degrees (6.1% variance). 

The mean and variance of the RMS pointing error of each axis is computed using LCP. 
The PDF and CDF of the RMS of PRX is plotted in figure 31(a), and of PRZ in figure 31(b). 
The PRX axis is both the largest contributor to the RSS (contributing 4.91 as compared 
to 2.91 deg) as well as to the variance. This relationship is made even more clear in the 
following analysis. 

A useful way to determine the relative importance of uncertain parameters is to look 
at the contribution of each to the total variance. This quantity is plotted for the RMS X 
and Z motions in figure 32. The top plot shows the mean value and 3-sigma bounds of 
each component of the RSS on the Y axis, for each component plotted on the X axis. The 
PRX RMS motion is clearly the dominant source. The bottom plot shows the normalized 
contributions of each uncertain parameter variance (color coded by the legend) to each 
component of the RSS. This result shows that the uncertainty in the PRX RMS variance 
is due primarily to variance in the Primary Appendage X-axis rotational stiffness. Since 
the PRX RMS is the largest contributor to the RSS, the same parameter dominates the 
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Figure 30: RSS pointing error computed using LCP. Top, CDF, bottom PDF. Red lines 
assume nominal design achieves a margin of 5%. 


uncertainty in the RSS as well. In the Z component, the active strut and Lexan center 
moduli are dominant contributors. 
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(a) PRX (b) PRZ 

Figure 31: CDF and PDF of PRX and PRZ RMS computed using LCP. 
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5.4 Change of Variables Results 

The output RSS is next examined using the exact transformation from the distributions 
of Lexan end and center moduli, and the Primary Appendage X rotational stiffness. The 
resulting CDF and PDF of J is plotted in figure 33. The red line is the exact PDF computed 
using CoV. The green line is the LCP solution. The blue bars are the Monte Carlo check 
computed using 2000 samples with Gaussian distribution. The CoV result shows that 
the true distribution of the RSS “tails” very slightly off to the right (high RSS) side. This 
results in a slightly decreased probability of success (74% versus 79%). This would, however, 
potentially create more difficulty in achieving a higher probability of success (above 99%) 
since the mean might need to be lowered significantly. Alternatively, lowering the sensitivity 
of the RSS to the uncertainty in the Primary Appendage Base X stiffness would reduce the 
variance of the RSS and should be “cheaper”. 


CO 0-6- 

rr 

n 

OT x 0.4- 


CDF of output RSS 




Figure 33: Comparison of CoV (red), LCP (green), and Monte Carlo (blue bars) for 3 
significant parameters. 

As a final check, a Monte Carlo analysis under uncertainty in all 8 parameters is per- 
formed, and compared to the CoV solution (figure 34). The Monte Carlo results do not 
show any dramatic differences compared to the CoV solution, verifying that the identified 
critical parameters are, in fact, the primary contributors to variance of the metric. 

The relative computation times of the LCP, CoV, and Monte Carlo methods are of 
interest. The computation times for the 3-parameter solution are given in table 5. The LCP 
solution is the quickest method, and also provides useful information about the relative 
importance of different parameters. However, it cannot accommodate nonlinear and/or 
non-monotonic cost functions, or non- Gaussian parameter distributions. The Change of 
Variables procedure takes substantial computational effort (not only to compute but also 
to regularize, or map onto a uniform grid, the samples). However, it gives the exact shape 
of the PDF/CDF for nonlinear/non-monotonic functions and arbitrary input distributions, 
and is efficient to compute the tails of the distribution. The Monte Carlo solution time 
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Figure 34: Comparison of CoV solution for 3 significant parameters (red) with Monte Carlo 
for all 8 parameters (blue bars). 


will increase moderately as the number of parameters increases. However, the number of 
computations needed to accurately map the tails may be problematic. 


Table 5: Relative computation times 


3 Parameter Solution 

Solution Points 

Time [sec] 

Linear Covariance Propagation 

i 

2.76 sec 

Change of Variables 

1000 

285.2 sec 



(+1331.4 sec to regularize) 

Monte Carlo 

2000 

467 sec 


5.5 Deliverables 

The Uncertainty Toolbox contains several uncertainty propagation tools, plotting functions, 
and analysis functions. A Command Reference Manual (DOCS ..Uncertainty .Reference) is 
also posted at the DOCS Web server. 

5.6 MACE Example Conclusions 

• Uncertainty analysis was performed for a realistic sized problem given: 

— a flight vehicle with flexible structural dynamics 

— integrated structural /electronic dynamics 

— closed loop (servo controlled gimbals and reaction wheels) 
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Table 6: MATLAB Deliverables 



Function 

Description 

propagation 

CoV_lD 

1-parameter Change of Variables 


CoV_2D 

2-parameter Change of Variables 


CoV_3D 

3-parameter Change of Variables 


lcp 

Linear Covariance Propagation 


montecarlo 

Monte Carlo simulation 

analysis 

binjnontecarlo 

Create histogram 


pdf _Ps 

Find Probability of Success from a ID PDF 


pdf .statistics 

Find mean/variance from ID PDF 

plot 

plot_pdf 

Plot ID PDF 


plot-gaussian 

Plot a Gaussian distribution 


pi 0 1 -gaus s i an_cont r ibut i ons 

Plot contribution of each parameter to total variance 


pi 0 1 jnont e c ar 1 0 

Plot Monte Carlo histogram 


• Large parameter variance (19%) results in only moderate objective variance (6.1%) 

- Large parameter variance was required to match FRF of model to data 

- Mode frequencies don’t have large influence on RSS (energy is still there) 

- the cost is nearly linear 

• Stability is a function of modal frequencies. 

— Should be more sensitive to critical parameters 

5.7 Independence of Output PDF 

The change-of- variables solution approach is powerful, but cannot be applied to practical 
systems for more than three parameters at a time. The question arises whether the depen- 
dence of the cost function to each parameter could be nearly independent. This would allow 
PDFs from separate, computationally tractable analyses to be combined. For independent 
random variables ( 21 , 2 : 2 ) the PDF of the joint distribution is the product of the individual 
distributions, 


fx 1,X2 ( 21 , 22 ) — fx 1 (xi) • f X2 {x 2 ) 


( 22 ) 


Unfortunately test cases showed that the coupling of parameters through a function 
invalidated the independence property. 
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6 Conclusion 


The objective of this research was to develop design methods and tools to the aerospace 
vehicle design process which take into account lifecycle uncertainties. The MIT SSL met 
this objective by investigating approaches to uncertainty identification and propagation, by 
using component testing to provide quantitative uncertainty descriptions, and by testing an 
approach to uncertainty analysis on the MACE testbed. 

MACE was used as a example of a flexible, high performance vehicle. New models of 
the structure have been built in order to recreate the “preliminary” design stage. This 
allows uncertainty analyses to be performed on an early-stage model, with the hardware 
available for analysis validation. A listing of potential uncertainty sources was provided, 
and engineering insight was used to select those sources most likely to affect the output 
performance. 

The benefits of uncertainty testing became very apparent. Non-linearity in the strut- 
node interfaces was considered a very likely source of uncertainty before testing showed it 
was not present. Testing saved time and resources on attempting to model an uncertainty 
that wasn’t there. Likewise, parametric testing for modulus, rotary inertias and stiffness 
provided data-based distributions in place of educated guesses. 

The change-of- variables propagation routine was developed. Although only capable of 
propagating at most three inputs through the system, it provides greater detail at the tails of 
the output distributions than linear covariance propagation or Monte Carlo sampling could. 
The routines were coded into MATLAB and have expanded the capabilities of DOCS. 

Lastly, all of these research directions were combined in an uncertainty analysis of 
MACE. Several propagation routines were used to identify critical uncertainty sources, 
and propagate these through the model. This showed that such uncertainty analyses are 
possible for realistically-sized problems. 
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7 Mapping Research to the SOW 


Statement of Work 

1) Use the DOCS toolset to conduct design 
optimizations to not only maximize per- 
formance but augment that design process 
(using MIT’s sensitivity and uncertainty 
analyses) with an ability to also minimize 
sensitivity to design error or uncertainty 
(maximize performance while minimizing 
error bars). 


2) Use Markov reliability analysis tools to de- 
velop a methodology for reliability opti- 
mization. If improvements in mean time 
to failure as a function of funds spent, and 
the cost of redundancy in design, can be 
modeled, optimal control techniques can be 
used to identify where funding can be most 
effectively spent in order to maximize life- 
cycle productivity in the presence of fail- 
ures. 

3) Apply these tools to two different 
Aerospace examples. In the first, fa- 
tigue life of aerodynamic surfaces will 
be used as an example for maximizing 
fatigue life while minimizing sensitivity to 
uncertainty. The second example will be 
interferometry where optical metrics will 
be optimized. 


Work Performed 

The DOCS toolset was used in analyz- 
ing the MACE model, including model up- 
dating and the uncertainty analysis. The 
design process augmented with measured 
parametric input bounds, along with the 
Change-of- Variables MATLAB codes (new 
addition to DOCS) plus other supporting 
codes. Sensitivity to uncertainty is reduced 
by identifying and characterizing uncertain 
parameters, and through the MACE un- 
certainty analysis. Importance of testing 
shown by elimination of suspected uncer- 
tainties (strut-node non-linearity). 

Method for reliability optimization devel- 
oped and published in Wertzs Masters the- 
sis [21]. 


Tools applied to two different aerospace ex- 
amples. Reliability applied to a distributed 
node network in Wertz’s Masters thesis. In 
place of aerodynamic surfaces, the uncer- 
tainty analysis applied to the pointing er- 
ror (akin to optical pointing metrics) of a 
flexible controls platform (MACE). Using 
MACE allowed direct testing of hardware 
which was a great benefit in examining the 
uncertainties between model and data. 
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